function fig4_18

%  solve IVP using MATLAB routines
%         y' = f(t,y)  with   y(0) = y0
%  where y = (y1, y2, y3 , ..., yn) is an n-vector

clf
% get(gcf)
set(gcf,'Position', [8 293 733 1067])
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)
hold on

tmax=160;
N=200;

%  initial values
y0=zeros(4*N-2,1);
y0(2*N-1)=1;

%  calculate solution using a MATLAB routine
[t,y] = ode45(@rhs,[0 tmax],y0);  titleX='ode45';
%[t,y] = ode23s(@rhs,[0 tmax],y0);  titleX='ode23s';

nt=size(t);

for it=1:nt(1)
    if t(it)<20
        nt1=it ;
    end
    if t(it)<40
        nt2=it ;
    end
    if t(it)<80
        nt3=it ;
    end
end
nt4=nt(1);

t(nt1)
t(nt2)
t(nt3)
t(nt4)

subaxis(5,1,1,1,'MT',0.005,'MB',0.05,'MR',0.018,'ML',0.09,'P',0.004)

jj=0;
sol=zeros(2*N-1,1);
for j=1:2:4*N-2
    jj=jj+1;
    sol(jj)=y(1,j);
end

x=linspace(-N+1,N-1,2*N-1);
hold on
plot(x,sol,'--r')
plot(x,sol,'bo','MarkerSize',2,'LineWidth',1)
grid on
box on
ax = gca;
ax.XAxis.TickLabelInterpreter = 'latex';
ax.YAxis.TickLabelInterpreter = 'latex';
set(gca,'FontSize',13)
ylabel('$y$\,-\,axis','Interpreter','latex','FontSize',14)
%t0=round(10*t(nt1))/10;
t0=round(t(1));
say=['$t\,=\,\,$', num2str(t0)];
text(45,0.85,say,'Interpreter','latex','FontSize',14)
axis([-N N -0.01 1]);

subaxis(5,1,1,2);

jj=0;
sol=zeros(2*N-1,1);
for j=1:2:4*N-2
    jj=jj+1;
    sol(jj)=y(nt1,j);
end

x=linspace(-N+1,N-1,2*N-1);
hold on
box on
plot(x,sol,'--r')
plot(x,sol,'bo','MarkerSize',2,'LineWidth',1)
grid on
ylabel('$y$\,-\,axis','Interpreter','latex','FontSize',14)

set(gca,'ytick',[-0.1 0 0.1 0.2])
set(gca,'YTickLabel',{'$-0.1$';'0';'0.1';'$0.2$'})

ax = gca;
ax.XAxis.TickLabelInterpreter = 'latex';
ax.YAxis.TickLabelInterpreter = 'latex';
set(gca,'FontSize',13)

%t0=round(10*t(nt1))/10;
t0=round(t(nt1));
say=['$t\,=\,\,$', num2str(t0)];
text(45,0.15,say,'Interpreter','latex','FontSize',14)
axis([-N N -0.15 0.2]);

subaxis(5,1,1,3);

jj=0;
sol=zeros(2*N-1,1);
for j=1:2:4*N-2
    jj=jj+1;
    sol(jj)=y(nt2,j);
end

x=linspace(-N+1,N-1,2*N-1);
hold on
plot(x,sol,'--r')
plot(x,sol,'bo','MarkerSize',2,'LineWidth',1)
grid on
box on
set(gca,'ytick',[-0.1 0 0.1 0.2])
set(gca,'YTickLabel',{'$-0.1$';'0';'0.1';'$0.2$'})

ax = gca;
ax.XAxis.TickLabelInterpreter = 'latex';
ax.YAxis.TickLabelInterpreter = 'latex';
set(gca,'FontSize',13)
ylabel('$y$\,-\,axis','Interpreter','latex','FontSize',14)
t0=round(t(nt2));
say=['$t\,=\,\,$', num2str(t0)];
text(45,0.15,say,'Interpreter','latex','FontSize',14)
axis([-N N -0.15 0.2]);

subaxis(5,1,1,4);

jj=0;
sol=zeros(2*N-1,1);
for j=1:2:4*N-2
    jj=jj+1;
    sol(jj)=y(nt3,j);
end

x=linspace(-N+1,N-1,2*N-1);
hold on
plot(x,sol,'--r')
plot(x,sol,'bo','MarkerSize',2,'LineWidth',1)
grid on
box on
set(gca,'ytick',[-0.1 0 0.1 0.2])
set(gca,'YTickLabel',{'$-0.1$';'0';'0.1';'$0.2$'})

ax = gca;
ax.XAxis.TickLabelInterpreter = 'latex';
ax.YAxis.TickLabelInterpreter = 'latex';
set(gca,'ytick',[-0.1 0 0.1 0.2])
set(gca,'YTickLabel',{'$-0.1$';'0';'0.1';'$0.2$'})

ax = gca;
ax.XAxis.TickLabelInterpreter = 'latex';
ax.YAxis.TickLabelInterpreter = 'latex';
set(gca,'FontSize',13)
ylabel('$y$\,-\,axis','Interpreter','latex','FontSize',14)
t0=round(t(nt3));
say=['$t\,=\,\,$', num2str(t0)];
text(45,0.15,say,'Interpreter','latex','FontSize',14)
axis([-N N -0.15 0.2])


subaxis(5,1,1,5);

jj=0;
sol=zeros(2*N-1,1);
for j=1:2:4*N-2
    jj=jj+1;
    sol(jj)=y(nt4,j);
end

x=linspace(-N+1,N-1,2*N-1);
hold on
plot(x,sol,'--r')
plot(x,sol,'bo','MarkerSize',2,'LineWidth',1)
grid on
box on
t0=round(t(nt4));
say=['$t\,=\,\,$', num2str(t0)];
text(45,0.15,say,'Interpreter','latex','FontSize',14)
axis([-N N -0.15 0.2])

set(gca,'ytick',[-0.1 0 0.1 0.2])
set(gca,'YTickLabel',{'$-0.1$';'0';'0.1';'$0.2$'})

ax = gca;
ax.XAxis.TickLabelInterpreter = 'latex';
ax.YAxis.TickLabelInterpreter = 'latex';
set(gca,'FontSize',13)
ylabel('$y$\,-\,axis','Interpreter','latex','FontSize',14)
xlabel('$n$\,-\,axis','Interpreter','latex','FontSize',14)

%set(gca,'xtick',[0 100 200 300 400 500]);

grid on
box on
set(gca,'FontSize',14)
hold off

%exportgraphics(gcf,'/Users/mark/Desktop/chain2.pdf')

%  define f(t,y)
function dy=rhs(t,y)
% y(odd) = y_i	y(even) = v_i
N=200;
a=1;
jj=1;
dy=zeros(4*N-2,1);
dy(jj)=y(jj+1);
dy(jj+1)=a*(y(jj+2)-2*y(jj));
for j=2:2*N-2
    jj=jj+2;
    dy(jj)=y(jj+1);
    dy(jj+1)=a*(y(jj+2)-2*y(jj)+y(jj-2));
end
jj=jj+2;
dy(jj)=y(jj+1);
dy(jj+1)=a*(-2*y(jj)+y(jj-2));
%pause





















